Custom functions
We also provide some additional helper functions to calculate effect
sizes, process data, and visualise our results. The most straightforward
way to use these custom functions is to run the code chunk below.
Alternatively, paste the code into the console and hit
Enter to have R ‘learn’ these custom
functions.
If you want to use these custom functions in your own data, you will
need to change the variable names according to your own data (check out
the R code and you will see what we mean).
# calculate effect sizes ----
## skewness ----
calc.skewness <- function(x, output = "est") {
n <- length(x)
if (output == "est") { # skewness estimate
(sqrt(n * (n - 1)) / (n - 2)) *
(((1 / n) * sum((x - mean(x)) ^ 3)) /
(((1 / n) * sum((x - mean(x)) ^ 2)) ^ (3/2)))
} else if (output == "var") { # skewness sampling variance
(6 * n * (n - 1)) /
((n - 2) * (n + 1) * (n + 3))
}
}
## kurtosis ----
calc.kurtosis <- function(x, output = "est") {
n <- length(x)
if (output == "est") { # kurtosis estimate
((((n + 1) * n * (n - 1)) / ((n - 2) * (n - 3))) *
(sum((x - mean(x)) ^ 4) / (sum((x - mean(x)) ^ 2) ^ 2))) -
(3 * ((n - 1) ^ 2) / ((n - 2) * (n - 3)))
} else if (output == "var") { # kurtosis sampling variance
(24 * n * ((n - 1) ^ 2)) /
((n - 3) * (n - 2) * (n + 3) * (n + 5))
}
}
## Zr ----
r.to.zr <- # Zr estimate
function(r) {
0.5 * log((1 + r) / (1 - r))
}
zr.variance <- # Zr variance
function(n) {
1 / (n - 3)
}
## other effect sizes (lnRR and lnVR) ----
calc.effect <- function(data = raw_data,
m) { # calculates other already established effect size statistics
escalc(measure = m,
m1i = data$mean_male,
m2i = data$mean_female,
sd1i = data$sd_male,
sd2i = data$sd_female,
n1i = data$n_male,
n2i = data$n_female,
var.names = c(paste0(m,
"_est"),
paste0(m,
"_var")))
}
# processing functions ----
process.ind_effects <- function(chosen_trait = "fat_mass",
measure = "KU_delta") {
ind_effects <-
df_meta_analysed %>%
filter(trait_name == chosen_trait,
phenotyping_center %in% c("CCP-IMG",
"HMGU",
"JAX",
"MRC H",
"TCP")) %>%
mutate(type = "individual") %>%
select(phenotyping_center,
strain_fig,
n = n_total,
est = paste0(measure, "_", "est"),
var = paste0(measure, "_", "var"),
lower = paste0(measure, "_", "lower"),
upper = paste0(measure, "_", "upper"))
model <- rma.mv(data = ind_effects,
yi = est,
V = var,
test = "t",
random = list(~ 1|phenotyping_center,
~ 1|strain_fig))
df_model <- data.frame(trait_name = chosen_trait,
est = model$beta[1],
var = model$se ^ 2,
lower = model$ci.lb,
upper = model$ci.ub,
phenotyping_center = "Mean",
strain_fig = "ES")
ind_effects %>%
bind_rows(df_model) %>%
mutate(est_type = measure,
centre_and_strain = factor(paste0(phenotyping_center,
"\n",
strain_fig))) %>%
mutate(centre_and_strain = factor(centre_and_strain,
levels = c("Mean\nES",
rev(levels(centre_and_strain)[-6]))))
}
process.cor_effects <- function(chosen_trait_1 = "fat_mass",
chosen_trait_2 = "heart_weight") {
df_effects_cor <-
df_raw %>%
filter(trait_name %in% c(chosen_trait_1,
chosen_trait_2),
phenotyping_center %in% c("CCP-IMG",
"HMGU",
"JAX",
"MRC H",
"TCP")) %>%
pivot_wider(id_cols = c(specimen_id,
strain_fig,
phenotyping_center,
sex),
names_from = trait_name) %>%
clean_names() %>%
drop_na() %>%
group_by(strain_fig,
phenotyping_center,
sex) %>%
group_modify(~ correlate(.x)) %>%
drop_na(all_of(chosen_trait_2)) %>%
ungroup() %>%
left_join(df_raw %>%
filter(trait_name %in% c(chosen_trait_1,
chosen_trait_2),
phenotyping_center %in% c("CCP-IMG",
"HMGU",
"JAX",
"MRC H",
"TCP")) %>%
pivot_wider(id_cols = c(specimen_id,
strain_fig,
phenotyping_center,
sex),
names_from = trait_name) %>%
clean_names() %>%
drop_na() %>%
group_by(strain_fig,
phenotyping_center,
sex) %>%
summarise(n = n())) %>%
rename(r_est = chosen_trait_2) %>%
mutate(zr_est = r.to.zr(r_est),
zr_var = zr.variance(n)) %>%
select(- c(4:6)) %>%
pivot_wider(names_from = sex,
values_from = c(n,
zr_est,
zr_var)) %>%
mutate(delta_zr_est = zr_est_male - zr_est_female,
delta_zr_var = zr_var_male + zr_var_female,
delta_zr_upper = delta_zr_est +
qt(0.975, n_male + n_female - 2) *
sqrt(delta_zr_var),
delta_zr_lower = delta_zr_est -
qt(0.975, n_male + n_female - 2) *
sqrt(delta_zr_var))
mlma_zr <-
rma.mv(data = df_effects_cor,
yi = delta_zr_est,
V = delta_zr_var,
test = "t",
random = list(~ 1|phenotyping_center,
~ 1|strain_fig))
df_model <- data.frame(delta_zr_est = mlma_zr$beta[1],
delta_zr_lower = mlma_zr$ci.lb,
delta_zr_upper = mlma_zr$ci.ub,
phenotyping_center = "Mean",
strain_fig = "ES")
df_effects_cor %>%
bind_rows(df_model) %>%
mutate(centre_and_strain = factor(paste0(phenotyping_center,
"\n",
strain_fig))) %>%
mutate(centre_and_strain = factor(centre_and_strain,
levels = c("Mean\nES",
rev(levels(centre_and_strain)[-5]))))
}
# visualisation functions ----
caterpillar.custom <-
function(chosen_trait = "fat_mass",
measure = "KU_delta") {
plot <-
process.ind_effects(chosen_trait = chosen_trait,
measure = measure) %>%
ggplot(aes(y = centre_and_strain,
x = est,
xmax = upper,
xmin = lower,
shape = strain_fig,
col = phenotyping_center)) +
geom_pointrange() +
geom_vline(xintercept = 0,
linetype = "dotted") +
theme_classic() +
theme(legend.position = "none",
axis.text.y = element_blank(),
axis.title.y = element_blank(),
plot.tag.position = c(0.15, 0.98))
if (measure == "ROM") {
plot +
labs(x = "lnRR") +
scale_x_continuous(limits = c(-0.51, 0.51),
breaks = c(-0.5, 0, 0.5)) +
theme(axis.title.x = ggtext::element_markdown(face = "italic"))
} else if (measure == "VR") {
plot +
labs(x = "lnVR") +
scale_x_continuous(limits = c(-1, 1),
breaks = c(-1, 0, 1)) +
theme(axis.title.x = ggtext::element_markdown(face = "italic"))
} else if (measure == "SK_delta") {
plot +
labs(x = "Δ*sk*") +
scale_x_continuous(limits = c(-2.1, 2.1),
breaks = c(-2, 0, 2)) +
theme(axis.title.x = ggtext::element_markdown())
} else if (measure == "KU_delta") {
plot +
labs(x = "Δ*ku*") +
scale_x_continuous(limits = c(-15, 15),
breaks = c(-15, 0, 15)) +
theme(axis.title.x = ggtext::element_markdown())
}
}
ridgeline.custom <- function(chosen_trait = "fat_mass") {
df_raw %>%
filter(trait_name == chosen_trait,
phenotyping_center %in% c("CCP-IMG",
"HMGU",
"JAX",
"MRC H",
"TCP")) %>%
add_row(phenotyping_center = "Mean",
strain_fig = "ES") %>%
mutate(centre_and_strain = factor(paste0(phenotyping_center,
"\n",
strain_fig))) %>%
mutate(centre_and_strain = factor(centre_and_strain,
levels = c("Mean\nES",
rev(levels(centre_and_strain)[-5]))),
value_s = scale(value)) %>%
ggplot(aes(x = value_s,
y = centre_and_strain,
fill = sex,
linetype = sex)) +
stat_slab(scale = 0.7,
alpha = 0.4,
linewidth = 0.6,
col = "black") +
scale_fill_manual(values = c("white",
"black")) +
scale_linetype_manual(values = c("solid",
"dashed")) +
labs(x = paste0(str_to_sentence(str_replace_all(chosen_trait,
"_",
" ")),
"\n(scaled)"),
y = "Phenotyping centre and mice strain") +
theme_classic() +
theme(legend.position = "none",
axis.title.x = element_text(size = 12,
margin = margin(t = 0.2,
unit = "cm")),
axis.title.y = element_text(size = 12,
margin = margin(r = 0.2,
unit = "cm")),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10),
plot.tag.position = c(0.53, 0.98))
}
cor.caterpillar.custom <-
function(chosen_trait_1 = "fat_mass",
chosen_trait_2 = "heart_weight") {
process.cor_effects(chosen_trait_1 = chosen_trait_1,
chosen_trait_2 = chosen_trait_2) %>%
ggplot(aes(y = centre_and_strain,
x = delta_zr_est,
xmax = delta_zr_upper,
xmin = delta_zr_lower,
shape = strain_fig,
col = phenotyping_center)) +
geom_pointrange() +
geom_vline(xintercept = 0,
linetype = "dotted") +
labs(y = "Phenotyping centre and mice strain",
x = "Δ*Zr*",
shape = "Strain") +
scale_x_continuous(limits = c(-1, 1),
breaks = c(-1, 0, 1)) +
theme_classic() +
theme(legend.position = "none",
axis.title.x = ggtext::element_markdown(size = 12,
margin = margin(t = 0.2,
unit = "cm")),
axis.title.y = element_text(size = 12,
margin = margin(r = - 0.1,
unit = "cm")),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10),
plot.tag.position = c(0.3, 0.99))
}
cor.plot.custom <-
function(chosen_trait_1 = "fat_mass",
chosen_trait_2 = "heart_weight",
chosen_lims = c(-3, 5)) {
df_cor <-
df_raw %>%
filter(trait_name %in% c(chosen_trait_1,
chosen_trait_2),
phenotyping_center %in% c("CCP-IMG",
"HMGU",
"JAX",
"MRC H",
"TCP")) %>%
pivot_wider(id_cols = c(specimen_id,
strain_fig,
phenotyping_center,
sex),
names_from = trait_name) %>%
clean_names() %>%
drop_na() %>%
mutate(centre_and_strain = factor(paste0(phenotyping_center,
strain_fig))) %>%
mutate(centre_and_strain = factor(centre_and_strain,
levels = rev(levels(centre_and_strain))),
trait_1_s = scale(get(chosen_trait_1))[,1],
trait_2_s = scale(get(chosen_trait_2))[,1])
plot_list <- list()
for (i in 1:length(levels(df_cor$centre_and_strain))) {
level_i <- sort(levels(df_cor$centre_and_strain))[i]
plot <-
df_cor %>%
filter(centre_and_strain == level_i) %>%
ggplot(aes(x = trait_1_s,
y = trait_2_s,
shape = sex,
linetype = sex)) +
geom_point(
alpha = 0.008,
) +
geom_abline(intercept = 0,
slope = 1,
linewidth = 0.5,
linetype = "dotted") +
geom_smooth(method = "lm",
se = F,
col = "black") +
scale_shape_manual(values = c(3, 4)) +
scale_linetype_manual(values = c("solid",
"dashed")) +
scale_x_continuous(limits = chosen_lims) +
scale_y_continuous(limits = chosen_lims) +
labs(x = paste0(str_to_sentence(str_replace_all(chosen_trait_1,
"_",
" ")),
"\n(scaled)"),
y = paste0(str_to_sentence(str_replace_all(chosen_trait_2,
"_",
" ")),
" (scaled)")) +
theme_classic() +
theme(legend.position = "none",
plot.tag.position = c(0.05, 0.91),
axis.title.x = element_text(size = 12,
margin = margin(t = 0.2,
unit = "cm")),
axis.title.y = element_text(size = 12,
margin = margin(r = 0.2,
unit = "cm")),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10))
if (i != 6) {
plot <-
plot +
theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.line.x = element_blank(),
axis.ticks.x = element_blank())
}
plot_list[[i]] <- plot
}
return(plot_list)
}